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A  three-dimensional  numerical  solver  is  developed  to  model  complex  transport  processes  inside  all  com¬ 
ponents  of  a  solid  oxide  fuel  cell  (SOFC).  An  initial  assessment  of  the  accuracy  of  the  model  is  made  by 
comparing  a  numerically  generated  polarization  curve  with  experimental  results.  Sensitivity  derivatives 
of  objective  functions  representing  the  cell  voltage  and  the  concentration  polarization  are  obtained  with 
respect  to  the  material  properties  of  the  anode  and  the  cathode  using  discrete  adjoint  method.  Imple¬ 
mentation  of  the  discrete  adjoint  method  is  validated  by  comparing  sensitivity  derivatives  obtained  using 
the  adjoint  technique  with  results  obtained  using  direct-differentiation  and  finite-difference  methods. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Research  and  development  of  alternative  energy  producing 
devices  has  been  gaining  momentum  in  recent  years  and  solid  oxide 
fuel  cells  (SOFCs)  are  considered  as  one  of  the  most  promising 
emerging  technologies.  Experimental  and  numerical  approaches 
have  been  undertaken  by  several  researchers  [1-8]  to  study  the 
behavior  of  SOFC  Some  of  the  advantages  of  the  numerical  sim¬ 
ulations  over  the  experiments  are  the  cost  effectiveness  and  the 
fact  that  the  simulations  provide  a  wealth  of  data  that  is  difficult 
or  impossible  to  obtain  experimentally  and  can  be  used  to  per¬ 
form  in-depth  analysis  of  the  SOFC  unit/ system.  Although  SOFC  are 
still  in  the  developmental  stage,  numerical  simulations  can  con¬ 
tribute  greatly  toward  better  designs  that  can  produce  more  power, 
increased  efficiency  and  extended  life  expectancy  of  various  SOFC 
components. 

To  date,  numerical  simulations  have  been  primarily  focused 
on  analysis  of  fuel  cells  or  fuel  cell  components,  without  strong 
emphasis  on  utilizing  the  simulations  in  a  design  optimization 
environment.  In  particular,  optimization  procedures  that  may  be 
efficiently  used  for  a  large  number  of  design  variables  have  not 
been  developed.  Because  of  the  emphasis  on  analysis  instead  of 
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design,  sensitivity  information  to  determine  the  effects  of  vari¬ 
ations  in  design  parameters  on  performance  has  been  primarily 
implemented  by  simply  changing  the  parameter  of  interest,  re¬ 
running  the  simulation,  and  comparing  the  results  with  those  from 
the  original  simulation  [1,5, 6, 8].  While  this  approach  can  be  used 
to  determine  the  effects  of  parameter  variations  on  fuel  cell  perfor¬ 
mance,  a  more  rigorous  approach  toward  optimization  would  likely 
lead  to  better  designs,  and  can  also  provide  improved  insight  into 
the  parameters  affecting  the  performance  of  the  fuel  cell.  For  SOFC 
problems,  example  cost  functions  that  can  be  used  for  improving 
performance  include  minimizing  temperature  variations,  obtaining 
equal  distribution  of  fuel  in  each  of  the  channels,  or  maximizing 
power.  Design  variables  may  be  related  to  the  shape/size  of  the  fuel 
channels,  electrodes,  electrolyte,  and  interconnect,  but  may  also 
be  coupled  to  the  stoichiometric  composition  of  fuel  or  material 
properties  such  as  the  porosity  or  tortuosity  of  the  electrodes. 

In  references  [9]  and  [10],  optimization  algorithms  have  been 
used  to  improve  the  performance  of  a  polymer-electrolyte- 
membrane  (PEM)  fuel  cell  using  four  design  variables,  where  the 
sensitivity  derivatives  used  for  the  optimization  algorithm  have 
been  obtained  using  a  finite-difference  approach.  While  finite 
differences  are  often  a  viable  means  for  computing  sensitivity 
derivatives,  this  method  can  be  computationally  restrictive  when 
a  sufficiently  large  number  of  design  variables  are  present.  In 
addition,  accurate  derivatives  can  sometimes  be  difficult  to  obtain 
using  finite  differences  because  of  subtractive  cancellation  errors 
[11  ],  which  occur  when  the  function  evaluations  in  the  numerator 
become  computationally  indistinguishable  [12]  when  very  small 
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Nomenclature 

B  permeability  (m2) 

e  internal  energy  (J  m-3 ) 

/  cost  function  (cost  function  depended) 

H  enthalpy  (J  kg-1) 

z  current  density  (A  m-2 ) 

z'o  exchange  current  density  (A  m-2 ) 

J  mass  flux  vector  (kg  nrr2  s_1 ) 

L  augmented  cost  function  (cost  function  dependent) 

m  mass  flow  rate  (kg  s_1 ) 

M  molecular  weight  (kg  kmol-1 ) 

ns  number  of  species 

P  pressure  (Nm-2) 

q  heat  flux  (J  m-2  s_1 ) 

Q.  solution  vector  (solution  variable  dependent) 

T  temperature  (K) 

u  x-velocity  component  (m  s_1 ) 

v  y-velocity  component  (m  s_1 ) 

w  z-velocity  component  (m  s_1 ) 

x,  y,  z  co-ordinate  system 
Xj  mole  fraction  of  zth  species 

Yt  mass  fraction  of  zth  species 

Greek  symbols 

V  control  volume  (m3) 

design  variable  vector  (design  variable  dependent) 
X  grid  vector  (m) 

s  porosity 

0  electric  potential  (V) 

r)  activation  polarization  (V) 

k  tortuosity 

A  costate  variable  vector  (cost  function  dependent) 

/x  molecular  viscosity  (kg ms-1) 

p  mass  concentration  (kg  nrr3 ) 

r  viscous  flux  (kg  m_1  s-2 ) 

Constants 

F  Faraday  constant  96484.56  (A  s  mol-1 ) 

Ru  universal  gas  constant  8314.4  (J  kmol-1  K-1 ) 

Indices 

a  anode 

c  cathode 

eff  effective 

z,  j  chemical  species 


perturbations  are  used.  By  using  a  discrete  adjoint  method,  sen¬ 
sitivity  derivatives  that  are  consistent  with  the  flow  solver  may  be 
obtained  for  use  in  a  design  optimization  environment.  A  particular 
strength  of  adjoint  methods  is  that  sensitivity  information  can  be 
obtained  with  a  computational  cost  that  is  only  weakly  dependent 
on  the  number  of  design  variables,  and  is  therefore  enabling  tech¬ 
nology  for  design  studies  where  many  design  variables  are  required. 

Secanell  et  al.  [13]  performed  gradient-based  optimization 
of  a  planar  self-breathing  polymer-electrolyte-membrane  fuel 
cell  cathode  using  a  two-dimensional  fuel  cell  model.  Sensi¬ 
tivity  derivatives  of  the  current  density  were  obtained  using  a 
direct-differentiation  method  with  respect  to  the  electrode  design 
parameters.  Although  direct  differentiation  is  an  accurate  method 
of  obtaining  sensitivity  derivatives,  application  of  this  method  for 
a  practical  design  problem  with  many  design  variables  is  computa¬ 
tionally  expensive. 


In  recent  years,  adjoint  methods  have  been  developed  and  uti¬ 
lized  for  numerical  simulations  in  the  aerodynamic  community 
for  sensitivity  analysis,  error  estimation,  and  adaptive  meshing 
[14-25].  Recently,  Kapadia  et  al.  [26]  performed  a  sensitivity  anal¬ 
ysis  for  a  three-dimensional  fuel  cell  type  geometry  where  the  cost 
function  was  based  on  the  requirement  of  equally  distributing  fluid 
through  the  channels.  This  numerical  experiment  included  all  the 
mesh  points  defining  the  surface  of  the  geometry  as  parameters, 
totaling  more  than  180,000  design  variables.  The  adjoint  method  is 
particularly  suited  to  this  class  of  problem  because  the  sensitivity 
derivatives  can  be  obtained  for  all  design  variables  with  the  compu¬ 
tational  cost  of  a  single  solution  of  the  non-linear  system  used  for 
analysis,  a  single  solution  of  the  linear  adjoint  system,  and  a  matrix- 
vector  multiply.  Although  the  numerical  model  used  in  this  problem 
did  not  include  diffusion  or  chemistry,  it  demonstrated  the  applica¬ 
bility  of  the  adjoint  method  for  solving  problems  with  many  design 
variables.  To  further  demonstrate  the  use  of  the  discrete  adjoint 
method  for  SOFC  applications,  Kapadia  et  al.  [26,27]  implemented 
the  adjoint  method  for  one-dimensional  [26]  and  two-dimensional 
[27]  SOFC  models.  In  these  studies,  Kapadia  et  al.  [26,27]  computed 
sensitivity  derivatives  of  several  cost  functions  reflecting  the  perfor¬ 
mance  of  SOFC  with  respect  to  geometric  and  material  properties 
of  the  fuel  cell. 

The  primary  goal  of  this  paper  is  to  formulate  and  develop 
adjoint  methodology  that  accounts  for  high-fidelity  physical  mod¬ 
eling  and  can  be  applied  to  practical  SOFC  design  applications 
in  three-dimensions.  For  demonstration  purposes,  the  geometry 
described  by  Wang  et  al.  [28]  has  been  chosen  due  to  the  availabil¬ 
ity  of  various  transport  coefficients  and  experimental  polarization 
curve.  To  verify  consistency  between  the  adjoint-based  sensitiv¬ 
ity  derivatives  and  the  flow  solver,  sensitivity  derivatives  obtained 
using  the  discrete  adjoint  method  are  compared  with  sensitivity 
derivatives  computed  using  the  direct-differentiation  and  finite- 
difference  methods. 

2.  Sensitivity  analysis 

2.1.  Discrete  adjoint  method 

The  goal  of  an  adjoint  method  is  to  determine  sensitivity  deriva¬ 
tives  that  can  be  used  in  a  formal  optimization  procedure  for 
minimizing  a  specified  cost  function,  which  is  indicative  of  the  per¬ 
formance  of  the  system.  A  general  optimization  procedure  begins  by 
first  defining  a  meaningful  cost  function  and  a  desired  set  of  design 
variables.  A  numerical  analysis  of  the  baseline  system  is  then  per¬ 
formed.  The  results  of  the  analysis  include  the  solution  variables  Q. 
of  the  discretized  partial  differential  equations,  which  are  subse¬ 
quently  used  to  determine  the  initial  cost.  Because  the  numerical 
analysis  involves  discretization  of  the  partial  differential  equations 
on  a  computational  mesh,  it  should  be  noted  that  Q.  represents  the 
vector  of  solution  variables  where  each  element  of  the  vector  is  rep¬ 
resentative  of  one  or  more  physical  variables  located  at  each  mesh 
point,  x- 

The  cost  function  may  have  an  explicit  dependence  on  the  vector 
of  design  variables,  but  will  also  have  an  implicit  dependence 
because  Qand  x  may  also  depend  on  the  design  variables.  Therefore, 
the  cost  function  is  typically  written  to  indicate  the  implicit  and 
explicit  dependence  on  the  design  variables  as, 

f=noiP),x{P),P)  (i) 

If  R  represents  the  vector  of  discrete  residuals  at  each  mesh 
point,  an  augmented  cost  function  I  can  be  defined  in  the  terms 
of  the  original  cost  function  and  the  vector  of  discrete  residuals  as 
following. 

U.QiP),  XlP),  P,  A)  =nQiP),  XlP),  P)  +  AtR(Q.(PI  XkP),  P)  (2) 
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In  Eq.  (2),  A  is  the  vector  of  Lagrange  multipliers  (also  known 
as  costate  variables).  Note  that  the  augmented  cost  function,  L,  is  a 
scalar  quantity  that  is  identical  to  the  original  cost  function/,  when 
K(Q)  is  zero,  indicating  that  the  steady-state  solution  is  obtained. 
Differentiating  the  augmented  cost  function  with  respect  to  each 
of  the  design  variables  yields  the  following  set  of  equations  for 
dljdp,  which  is  a  column  vector  where  each  element  represents 
the  derivative  of  the  augmented  cost  function  with  respect  to  a 
particular  design  variable. 


dL 

dp 


1 

dR 

T 

+ 

d/ 

T 

~dR~ 

\ 

.¥ 

_dP 

_9x_ 

A 


0) 


Because  the  elements  of  A  are  arbitrary,  the  second  term,  which 
involves  the  derivatives  of  the  dependent  variables  with  respect  to 
the  design  variables,  can  be  eliminated  by  solving  a  linear  system 
of  equations  for  the  costate  variables,  also  known  as  the  adjoint 
equation. 


Now,  R{QlP),xiP),P)  =  0 


(7) 


dR 

dR 

dQ 

dR 

dX 

dp 

+ 

dQ 

dp 

+ 

dx 

dp 

dR 

dQ 

dR 

" dR " 

dx 

dQ 

dp 

dp 

dx 

dp 

As  seen,  computation  of  dQJdp  is  an  essential  part  of  this 
method,  which  requires  the  solution  of  a  linear  system  of 
equations  for  each  design  variable.  This  requirement  makes  direct- 
differentiation  methods  computationally  expensive  for  problems 
with  many  design  variables.  However,  because  derivatives  of 
dependent  variables  with  respect  to  the  design  variables  are  com¬ 
puted  at  each  node  in  the  flowfield,  this  method  is  particularly 
useful  when  there  are  many  flowfield  constraints. 

Although  the  emphasis  of  this  paper  is  in  the  development  and 
use  of  an  adjoint-based  method,  the  direct-differentiation  approach 
has  also  been  used  and  results  will  be  presented  using  both  tech¬ 
niques. 


3.  Governing  equations  and  numerical  details 


dR 

9Q 


(4) 


Once  the  costate  variables  are  obtained,  the  derivatives  of  the 
cost  function  with  respect  to  all  the  design  variables  are  obtained 
using  a  matrix-vector  multiplication. 


dL 

dp 


dp\  dx 


+ 


'dR' 

T 

dX 

T 

~dR~ 

dp_ 

+ 

dP_ 

_9x_ 

A  (5) 


In  numerical  simulations,  the  largest  computational  cost  of  com¬ 
puting  sensitivity  derivatives  using  the  adjoint  equations  is  due  to 
the  solution  of  the  analysis  equations  and  the  adjoint  equation,  both 
of  which  are  independent  of  the  number  of  design  variables.  The 
only  dependency  on  the  number  of  design  variables  is  in  the  eval¬ 
uation  of  Eq.  (5),  which  is  generally  much  cheaper  to  compute  than 
either  the  analysis  or  adjoint  solutions. 

Note  that  the  terms  in  Eqs.  (2)-(5)  involve  differentiation  of  the 
discrete  residual  R ,  the  cost  function/,  and  the  computational  mesh 
X  with  respect  to  the  dependent  variables  Q,  the  design  variables  p, 
and  the  location  of  the  mesh  points  x •  Correct  implementation  of 
this  procedure  can  be  extremely  tedious  to  accomplish  by  hand  and 
the  resulting  code  can  be  difficult  to  maintain.  To  overcome  the  dif¬ 
ficulties  associated  with  hand  differentiation,  the  complex-variable 
technique  of  Burdyshaw  and  Anderson  [16]  and  Nielsen  and  Kleb 
[24]  has  been  used  for  evaluating  all  the  terms  in  the  matrices 
required  for  solving  the  adjoint  equations  and  for  evaluating  Eq.  (5) 
once  the  costate  variables  have  been  obtained.  Step-by-step  deriva¬ 
tion  of  complex-variable  technique  is  demonstrated  by  Kapadia  et 
al.  [26]  along  with  the  detailed  discussion  on  relative  benefits  and 
drawbacks  of  complex-variable  method  with  respect  to  automatic 
differentiation  [29,30]  and  finite-difference  methods. 


2.2.  Direct  differentiation 


As  described  earlier,  sensitivity  derivatives  can  also  be  computed 
using  the  direct-differentiation  method.  Derivation  of  this  method 
using  the  chain  rule  is  shown  in  Eqs.  (6)-(9)  below. 


[  dfmPfxiPIP)  1  [9/1  .  [9/1 

\  dp  J  Idpj  +  XdQ.) 


dx 

dp 


(6) 


3.1.  Governing  equations 

The  three-dimensional  SOFC  model  developed  in  this  paper 
is  essentially  an  extension  of  the  two-dimensional  SOFC  model 
presented  in  previous  work  [27]  and  also  includes  the  addition 
of  an  electric  potential  equation  that  governs  the  distribution  of 
electric  potential  and  current  density  in  the  flowfield.  The  three- 
dimensional  model  accounts  for  all  components  of  the  SOFC, 
including  the  anode,  cathode,  electrolyte,  interconnects,  and  the 
fuel  and  air  channels.  Note  that  the  model  is  not  limited  to  any  par¬ 
ticular  type  of  SOFC,  i.e.  planar  as  well  as  tubular  type  SOFC  can 
be  simulated  using  this  model.  In  the  presented  work,  simulations 
have  been  performed  on  a  planar  type  porous-electrode  supported 
SOFC  as  explained  by  Wang  et  al.  [28]. 

Geometrical  details  of  a  porous-electrode  supported  SOFC  [28] 
are  given  in  Fig.  1 .  The  figure  shows  the  front  view  of  the  actual 
geometry  used  in  the  numerical  simulation.  As  the  name  suggests, 
fuel  and  air  channels  are  bored  through  the  anode  and  cath¬ 
ode,  respectively.  A  very  thin  electrolyte  (0.05  mm)  is  sandwiched 
between  the  porous  electrodes  and  interconnects  are  located  at  the 
top  and  the  bottom  surfaces  of  the  SOFC  unit  in  Fig.  1.  The  geometry 
depicted  in  Fig.  1  is  extruded  in  the  z-direction  to  form  the  three- 
dimensional  geometry.  Dimensions  of  the  geometry  are  obtained 
from  literature  [28]  and  also  tabulated  in  Table  1. 

Fuel  mixture  and  air  enter  through  the  inlet  of  the  channels 
and  diffuse  inside  the  anode  and  cathode,  respectively.  Oxygen 
atoms  combine  with  the  electrons  and  get  converted  into  oxygen 
ions  (0.502  +  2e_  O2-).  These  oxygen  ions  migrate  through 

the  solid  electrolyte  to  reach  the  anode-electrolyte  interface 


Air  Interconnect  Cathode 
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Table  1 

Dimensions  of  porous-electrode  supported  SOFC  [28]. 


Length  (z-direction) 

60  mm 

Width  (x-direction) 

60  mm 

Anode  thickness  (y-direction) 

10  mm 

Cathode  thickness  (y-direction) 

10  mm 

Electrolyte  thickness  (y-direction) 

0.05  mm 

Interconnect  thickness  (y-direction) 

0.5  mm 

Channel  diameter  (x-direction) 

5.0  mm 

Distance  between  consecutive  channels  (x-direction) 

5.0  mm 

and  combine  with  hydrogen  to  generate  steam  and  release  elec¬ 
trons  (H2  +  O2-  H20  +  2e-).  Both  these  reactions  take  place 

inside  extremely  thin  layers  near  the  cathode-electrolyte  and  the 
anode-electrolyte  interfaces,  respectively.  Because  modeling  the 
details  of  the  interface  region  is  impractical  due  to  the  small  size, 
the  cumulative  effect  of  the  electrochemical  reactions  is  modeled 
as  a  jump  in  the  electric  potential. 

Governing  equations  for  the  mass,  momentum  and  energy 
conservation  are  solved  simultaneously  with  the  equation  govern¬ 
ing  the  electric  potential  in  the  numerical  model.  The  system  of 
equations  utilized  in  the  model  is  given  by  Eqs.  ( 10)-(  1 5 ),  which  rep¬ 
resent  the  conservation  statements  for  the  species  concentrations, 
momentum  (x,  y  and  z),  energy  and  current,  respectively. 


^|^  +  v.(eAV)  +  v.(7f)  =  si 

(10) 

(U) 

*^+V.(ep^).-«J+V^er,,-s!7 

(12) 

(13) 

^^+V.(eyoHV)  +  V.(^lHi) 

ns 

=  V  •  (efi rt)-  V  •  qeff  +  W<p  ■  [aVip) 

(14) 

V  ■  (ctV0)  =  0 

(15) 

Table  2 

Material  properties  of  various  components  of  SOFC  [28,31  ]. 

Electric  resistivity  of  anode  (£2m) 

2.98  xl0“5  exp(- 1392/7) 

Electric  resistivity  of  cathode  (£2m) 

8.11  x  10-5exp( -600/7) 

Electric  resistivity  of  interconnect  (ftm) 

6.41  x  10-8 

Ionic  resistivity  of  electrolyte  (ftm) 

2.94  x  10-5  exp(  10350/7) 

Thermal  conductivity  of  anode  (W  m_1 I<-1 ) 

6.23 

Thermal  conductivity  of  cathode  (W  m-1 I<-1 ) 

9.6 

Thermal  conductivity  of  interconnect  (W  m-1 K-1 ) 

9.6 

Thermal  conductivity  of  electrolyte  (W  m-1 K-1 ) 

2.7 

Porosity  of  anode 

0.38 

Porosity  of  cathode 

0.5 

Tortuosity  of  anode 

1.5 

Tortuosity  of  cathode 

1.5 

Permeability  of  anode  (m2 ) 

1.0  x  10-12 

Permeability  of  cathode  (m2) 

o 

X 

o 

Pore  diameter  of  anode  (m) 

2.0  x  10“6 

Pore  diameter  of  cathode  (m) 

2.0  x  10-6 

ing  activation  polarization,  concentration  polarization  and  ohmic 
polarization.  Noren  and  Hoffman  [32]  have  provided  extensive  dis¬ 
cussion  on  accurately  modeling  the  activation  polarization.  The 
SOFC  model  used  in  this  work  employs  the  Butler-Volmer  equation 
to  compute  activation  polarization  [28-29,32].  The  Butler-Volmer 
equation  can  be  written  as, 


i  =  io 


exp 


- exp  ((1  - 


(16) 


The  activation  polarization  is  denoted  by  r\act •  ol  is  the  charge 
transfer  coefficient  and  assumed  to  be  0.5  in  the  current  work.  ne 
represents  the  number  of  electrons  involved  in  the  electrochemical 
reaction,  which  is  2  in  the  current  simulation.  z0  is  the  exchange 
current  density  and  it  is  computed  using  Eqs.  (17)  and  (18)  for  the 
anode  and  cathode  [33],  respectively. 


Various  constants  in  the  above  equations  are  given  in  Table  3 
[33].  Once  the  values  of  a  and  ne  are  inserted  in  Eq.  (16),  the  activa¬ 
tion  polarization  can  be  computed  using  the  following  expression. 


Eqs.  (10)-(14)  are  modified  Navier-Stokes  equations  valid  for 
both  porous  and  fluid  regions.  Detailed  discussion  on  flux  formu¬ 
lation  for  these  equations  can  be  found  in  previous  work  [26,27]. 
Eq.  (15)  represents  the  electric  potential  equation.  Electric/ionic 
conductivity,  cr,  in  Eq.  (15)  is  a  strong  function  of  the  temperature. 
Expressions  describing  the  relationships  between  the  electric/ionic 
resistivity  (reciprocal  of  conductivity)  and  the  temperature  for  vari¬ 
ous  components  of  SOFC  are  presented  in  Table  2  [28,31  ]  along  with 
thermal  conductivities  and  other  material  properties  of  different 
components  of  the  SOFC. 

As  presented,  Eq.  (15)  is  an  elliptic  equation  contrary  to 
the  rest  of  the  governing  equations,  Eqs.  (10)-(14),  which  are 
hyperbolic-parabolic  equations.  Eq.  (15)  is  solved  in  the  entire 
domain  except  for  the  fuel  and  air  channels,  which  are  pure  fluid 
regions.  Boundary  conditions  utilized  while  solving  Eq.  (15)  are  the 
same  as  applied  by  Wang  et  al.  [28  ].  The  lower  boundary  of  the  SOFC 
shown  in  Fig.  1  (bottom  surface  of  the  interconnect)  is  assumed  to 
be  at  zero  potential  while  average  current  density  is  specified  on  the 
upper  boundary  of  the  interconnect  at  the  top  of  the  cell.  Thus,  the 
computed  potential  on  the  top  surface  of  the  interconnect  gives  the 
operating  voltage  of  the  SOFC.  Similar  boundary  conditions  have 
also  been  used  in  references  [1,6,28]. 

The  voltage  output  of  the  SOFC  strongly  depends  on  several 
irreversibilities  or  losses  encountered  in  the  flowfield  includ¬ 


Ohmic  polarization  is  a  direct  consequence  of  the  resistance 
offered  to  the  flow  of  electrons/ions  inside  various  components 
of  the  SOFC.  Voltage  drop  due  to  ohmic  resistance  is  directly  pro¬ 
portional  to  the  current  and  the  resistance.  The  effect  of  ohmic 
polarization  on  the  voltage  loss  is  directly  included  in  the  poten¬ 
tial  equation,  Eq.  (15),  through  the  electric  conductivity,  cr,  which  is 
the  reciprocal  of  the  electric  resistivity. 

Concentration  polarization  is  caused  by  reductions  in  the  con¬ 
centrations  of  the  reacting  species  at  the  interface  between  the 
electrodes  and  the  electrolyte.  The  effect  of  the  reduction  in  con¬ 
centrations  can  be  seen  from  the  well-known  Nernst  potential 
equation,  given  by  Eq.  (20).  Also,  exchange  current  densities  at  the 


Table  3 

Constants  used  to  compute  activation  polarization  [33]. 


a 

0.5 

ne 

2 

Za  (Am-2) 

5.5  x  108 

£c  (Am-2) 

7.0  x  108 

Eact,a  (Jkmol-1) 

1.0  x  108 

Eact,c  (J  kmol-1 ) 

1.2  x  108 

Pre/(Nm-2) 

101,325 
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anode-electrolyte  interface  and  the  cathode-electrolyte  interface, 
represented  by  Eqs.  (17)  and  (18),  respectively,  are  strongly  affected 
by  the  concentration  polarization. 

Eq.  (20)  computes  the  electromotive  force  (EMF)  or  electric 
potential  under  reversible  conditions,  i.e.  in  the  absence  of  acti¬ 
vation,  ohmic  or  any  other  losses. 

EMF  -  EMF°  +  ^  In  )  (20) 

2F  v  Ph2°  J 


Pi  = 


Pi 

Pref 


(21) 


The  electromotive  force  at  standard  pressure  is  given  by  EMF0. 
The  value  of  Pref  is  taken  as  one  atmosphere  in  above  equation. 

The  electrochemical  reaction  reduces  the  concentration  of  the 
reactants  and  increases  the  concentration  of  the  products  at  the 
electrode-electrolyte  interface.  Thus,  the  partial  pressures  of  the 
reactants  and  products  are  affected  in  the  same  manner.  This  will 
reduce  the  value  of  the  second  term  on  the  right-hand  side  of 
Eq.  (20)  thereby  affecting  the  EMF  of  the  cell  adversely.  Concen¬ 
tration  polarization  strongly  depends  on  the  material  properties 
of  the  electrodes  that  are  responsible  for  the  transport  (diffusion 
and  convection)  of  the  reactants  and  products,  to  and  from  the 
electrode-electrolyte  interface. 


^ applied 


3.2.  Boundary  and  interface  conditions 

Figs.  2  and  3  show  boundary  and  interface  conditions  applied  to 
the  computational  domain.  Fig.  3  shows  the  “YZ-plane”  extracted 
from  the  computational  geometry  shown  in  Fig.  1  along  the  stream- 
wise  direction.  As  presented,  the  plane  is  extracted  through  the 
mid-point  of  the  fuel  and  air  channels  such  that  all  components  of 
the  SOFC  are  visible.  Fig.  2  represents  the  front  view  of  the  compu- 


lapplied 


Fig.  2.  Boundary  and  interface  conditions  (front  view). 


tational  geometry  assuming  a  single  channel  instead  of  six  channels 
although  it  should  be  noted  that  all  six  channels  have  been  included 
in  the  analysis.  No-slip,  adiabatic  wall  boundary  conditions  are 
applied  at  the  top  wall,  bottom  wall  and  side  walls  as  shown  in 
Figs.  2  and  3.  As  mentioned  earlier,  fixed  potential  (0  =  0)  boundary 
condition  is  applied  at  the  bottom  wall.  The  top  wall  is  treated  by 
specifying  average  current  density  (i  =  iappiied )• 

Inflow  boundary  conditions  with  specified  mass  flow  rate  and 
species  mole  fractions  are  applied  at  both  fuel  and  air  channel  inlets. 
The  temperature  of  the  air  and  fuel  mixture  entering  at  their  respec¬ 
tive  channels  is  1273  K  [28].  Also,  both  channels  are  operating  at 
atmospheric  pressure.  Specified  back  pressure  outflow  conditions 
are  applied  at  both  air  and  fuel  channel  outlets. 

Several  transport  processes  take  place  at  the  anode-electrolyte 
and  the  cathode-electrolyte  interfaces  that  strongly  affect  the  over¬ 
all  behavior  of  the  SOFC.  The  conversion  of  oxygen  molecules  into 
oxygen  ions  at  the  cathode-electrolyte  interface  is  modeled  by 
applying  a  mass  flux  condition  at  the  interface  using  Faraday’s 
law.  In  Eq.  (22),  i  is  the  local  current  density  and  F  is  Faraday’s 
constant.  A  negative  sign  implies  that  the  flux  is  leaving  the 
interface. 

Jo2=-^M02  (22) 

Similarly,  the  following  mass  flux  conditions  for  hydrogen  and 
steam  are  applied  at  the  anode-electrolyte  interface. 

Jh2=-^Mh2  (23) 

Jh20  =  ^pMH2o  (24) 

3.3.  Solution  procedure 

Flowfield  variables  are  computed  using  an  unstructured, 

implicit,  finite-volume  solver.  The  solver  is  vertex  centered  and  the 
discrete  residual  at  each  node  is  computed  by  integrating  the  gov¬ 
erning  Eqs.  (10)-(15)  over  a  median  dual  control  volume.  Because 
a  steady-state  solution  is  the  primary  goal  of  the  current  work, 
time  accuracy  of  the  solution  is  sacrificed  by  allowing  local  time¬ 
stepping  to  accelerate  convergence. 
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Table  4 

Operating  conditions  utilized  in  polarization  curve. 
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*h2 

^h2o 

Xo2 

Xn2 

T(I<) 

P(  Nirr2) 

mjueii  kgs-1) 

mair(  kgs-1) 

0.9578 

0.0422 

0.198 

0.802 

1273 

101,325 

5.94  x  IQ"7 

2.15  x  10-5 

To  reduce  computer  time,  the  solution  is  obtained  using  multiple 
processors  utilizing  the  message  passing  interface  (MPI)  [34]  and 
necessary  grid  decomposition  is  achieved  using  METIS  [35].  Orig¬ 
inal  grids  are  generated  using  the  commercial  software  Gridgen 
[36]. 

An  implicit  Euler  scheme  is  used  to  solve  the  non-linear  sys¬ 
tem  as  given  by  Eqs.  (10)-(15).  A  flux-difference  splitting  scheme 
based  on  the  ROE  scheme  [37,38]  for  a  multi-component  mixture 
is  derived  to  model  the  convective  fluxes.  A  central-difference  for¬ 
mulation  is  used  to  compute  all  the  second-order  derivative  terms. 
Linear  systems  encountered  in  both  the  flowfield  and  sensitivity 
solvers  are  solved  using  the  GMRES  [39]  method. 

4.  Results  and  discussion 

4.1.  Polarization  curve  and  analysis  results 

As  mentioned  earlier,  the  numerical  simulations  presented  in 
this  paper  utilize  the  geometry  used  by  Wang  et  al.  [28].  One  of 
the  main  reasons  for  using  this  geometry  is  the  availability  of  an 
experimental  polarization  curve  [28]  that  can  be  used  to  validate 
the  three-dimensional  SOFC  model.  The  fuel  mixture  is  assumed  to 
contain  hydrogen  and  steam.  Air  is  modeled  as  a  mixture  of  oxygen 
and  nitrogen.  Species  mole  fractions  of  the  fuel  mixture  and  air 
entering  the  respective  channels  are  given  in  Table  4.  The  operating 
pressure,  temperature  and  mass  flow  rates  of  fuel  and  air  are  also 
given  in  Table  4. 

A  comparison  between  the  experimental  polarization  curve  [28] 
and  the  polarization  curve  obtained  using  the  numerical  model  is 
shown  in  Fig.  4.  As  can  be  seen,  the  overall  comparison  is  satis¬ 
factory.  The  numerical  tool  successfully  predicts  the  shape  of  the 
polarization  curve  and  obtains  results  that  are  within  two  percent 
of  the  experimental  data  at  low  current  densities  and  are  essentially 
identical  to  the  experimental  data  at  higher  current  densities.  As 


Fig.  4.  Polarization  curve. 


expected,  the  cell  voltage  reduces  with  increasing  current  density 
due  in  part  to  ohmic  losses  which  are  linearly  proportional  to  the 
current  density.  Also,  increases  in  current  draws  more  hydrogen  and 
oxygen  from  the  anode-electrolyte  and  the  cathode-electrolyte 
interfaces,  respectively  and  produces  more  steam.  This  reduces  the 
value  of  “EMF”  in  Eq.  (20),  i.e.  concentration  polarization  increases. 
Thus,  the  cumulative  effects  of  activation  polarization,  ohmic  polar¬ 
ization  and  concentration  polarization  are  evident  in  Fig.  4. 

Fig.  5  shows  the  surface  grid  used  in  the  current  simulation. 
The  figure  indicates  the  locations  of  the  planes  utilized  for  plot¬ 
ting  various  results.  Plane-A  passes  through  the  cathode  and  the 
air  channel  in  the  streamwise  direction.  Similarly,  Plane-B  passes 


Fig.  5.  Surface  grid. 
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Fig.  6.  Hydrogen  mole  fraction  inside  the  fuel  channel  and  the  anode. 
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Fig.  7.  Steam  mole  fraction  inside  the  fuel  channel  and  the  anode. 


through  the  anode  and  the  fuel  channel  in  the  streamwise  direction. 
Flow  directions  of  the  air  and  the  fuel  mixture  are  also  indicated  in 
Fig.  5.  Contours  of  hydrogen  and  steam  mole  fractions  are  plotted  on 
Plane-B  in  Figs.  6  and  7,  respectively.  Oxygen  mole  fraction  contours 
are  plotted  on  Plane-A  in  Fig.  8. 

Figs.  6-10  show  results  obtained  using  the  same  operating  con¬ 
ditions  as  used  in  the  polarization  curve.  Average  current  density  of 
4000  Am-2  is  applied  at  the  upper  boundary  of  the  interconnect. 
Note  that  the  current  density  is  computed  at  each  iteration  and 
it  varies  according  to  the  potential  gradient,  transport  properties 
and  operating  conditions  throughout  the  domain.  Fig.  6  shows  the 
contours  of  hydrogen  mole  fraction  plotted  in  the  Plane-B  shown 
in  Fig.  5.  Three  major  transport  processes  involving  hydrogen  take 
place  during  an  operation  cycle  of  SOFC:  (1 )  convection  of  hydrogen 
from  the  inlet  of  the  fuel  channel  to  the  outlet,  (2 )  diffusion  of  hydro¬ 
gen  from  the  fuel  channel  inside  the  anode  and  (3)  consumption 
of  hydrogen  due  to  the  electrochemical  reaction  occurring  at  the 
anode-electrolyte  interface.  The  combined  effects  of  these  trans¬ 
port  processes  are  evident  in  Fig.  6,  which  shows  a  reduction  of 
hydrogen  mole  fraction  along  the  flow  direction.  Careful  inspection 
of  Fig.  6  indicates  that  the  spot  with  the  lowest  mole  fraction  of 
hydrogen  is  located  at  the  upper  right  corner.  This  spot  is  located 
at  the  anode-electrolyte  interface  exactly  above  the  fuel  channel 
outlet. 

Fig.  7  shows  the  steam  mole  fraction  contours  plotted  in  the 
Plane-B  shown  in  Fig.  5.  The  top  boundary  in  Fig.  7  represents 
the  anode-electrolyte  interface.  As  described  earlier,  steam  is  pro¬ 
duced  due  to  the  electrochemical  reaction  at  the  anode-electrolyte 
interface  and  diffuses  inside  the  anode  and  eventually  into  the  fuel 
channel.  Thus,  the  mole  fraction  of  steam  increases  gradually  in 


the  streamwise  direction.  Mole  fraction  contours  of  steam  exhibit 
exactly  the  opposite  trend  as  shown  by  hydrogen  mole  fraction  con¬ 
tours  in  Fig.  6.  A  spot  with  the  maximum  steam  mole  fraction  is 
located  at  the  anode-electrolyte  interface,  just  above  the  fuel  chan¬ 
nel  outlet.  Nernst  potential,  represented  by  Eq.  (20)  reduces  along 
the  flow  direction  due  to  the  gradual  increase  in  steam  mole  fraction 
and  decrease  in  hydrogen  mole  fraction. 

Fig.  8  shows  the  oxygen  mole  fraction  contours  plotted  in  the 
Plane-A  shown  in  Fig.  5.  The  bottom  boundary  in  Fig.  8  represents 
the  cathode-electrolyte  interface.  Note  that  because  oxygen  acts  as 
a  reactant  in  the  electrochemical  reaction,  the  mole  fraction  of  oxy¬ 
gen  in  the  air  channel  reduces  along  the  flow  direction.  As  the  mass 
flow  rate  of  air  is  higher  than  the  mass  flow  rate  of  fuel,  the  over¬ 
all  reduction  in  oxygen  mole  fraction  is  less  than  the  reduction  in 
hydrogen  mole  fraction.  The  location  with  the  lowest  oxygen  mole 
fraction  is  the  cathode-electrolyte  interface  (lower  right  corner), 
just  below  the  air  channel  outlet. 

Fig.  9  shows  temperature  contours  plotted  on  the  outer  surfaces 
of  the  computational  domain.  Various  modes  of  energy  transfer 
that  contribute  to  the  temperature  distribution  include  convec¬ 
tion/diffusion/conduction  of  energy,  heat  generated  due  to  viscous 
stresses  and  most  importantly,  heat  generated  due  to  the  electro¬ 
chemical  reaction  at  the  anode-electrolyte  interface.  Also,  ohmic 
heating  contributes  to  the  rise  in  temperature.  Due  to  the  heat  gen¬ 
erated  during  the  electrochemical  reaction  at  the  anode-electrolyte 
interface,  there  is  a  gradual  increase  in  temperature  in  the  stream- 
wise  direction. 

The  extent  of  temperature  rise  depends  strongly  on  the  mass 
flow  rates  of  fuel  and  air.  Convective  cooling  increases  with  higher 
flow  rates  of  either  fuel  or  air  and  thus,  contributes  toward 
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Fig.  8.  Oxygen  mole  fraction  inside  the  air  channel  and  the  cathode. 


S.  Kapadia,  W.K.  Anderson  /  Journal  of  Power  Sources  189  (2009)  1074-1082 


1081 


1273  1285  1295  1305  1315  1325 
Fig.  9.  Temperature  contours  inside  the  computational  domain. 


Table  5 

Validation  of  adjoint  implementation  (Cost-1). 


D.V. 

Discrete  adjoint 

Direct  differentiation 

Finite  difference 

Sa 

-1.4136629036e-02 

-1.4136629036e-02 

-1.411798953e-02 

Ka 

-3.4924988954e-03 

-3.4924988954e-03 

-3.491000966e-03 

( r)a 

8.7578106870e+02 

8.7578106869e+02 

8.754119812e+02 

Sc 

2.7292696323e-03 

2.7292696322e-03 

2.761454211e-03 

KC 

-  1.4976041763e-03 

-  1.4976041763e-03 

-1.497321755e-03 

(r)c 

1.8945315028e+02 

1.8945315028e+02 

1.894163071  e+02 

Table  6 

Validation  of  adjoint  implementation  (Cost-2). 


D.V. 

Discrete  adjoint 

Direct  differentiation 

Finite  difference 

Sa 

-1.3870751906e-02 

-1.3870751907e-02 

-1.38668867e-02 

Ka 

-3.5959382450e-03 

-3.5959382450e-03 

-3.59483888e-03 

(r)a 

9.0264297426e+02 

9.0264297426e+02 

9.02379501  e+02 

Sc 

1.6938729757e-03 

1.6938729757e-03 

1.69193541e-03 

Kc 

8.7332085595e-05 

8.7332085595e-05 

8.75493837e-05 

(r)c 

-1.0164312515e+01 

-1.0164312515e+01 

-1.01922398e+01 

reduction  in  average  temperature  rise.  The  opposite  is  true  for  the 
lower  mass  flow  rates.  Current  density  is  another  major  contribu¬ 
tor  in  temperature  rise  due  to  the  fact  that  ohmic  heating  is  linearly 
proportional  to  the  current  density.  Also,  current  density  dictates 
the  extent  of  electrochemical  activity  at  the  electrode-electrolyte 
interface  that  is  responsible  for  generating  heat. 

Fig.  10  shows  current  density  vectors  plotted  in  the  inlet  plane 
shown  in  Fig.  1.  Note  that  large  variations  in  the  current  density  are 
observed  throughout  the  cell,  particularly  between  adjacent  chan¬ 
nels.  This  is  required  to  satisfy  the  current  conservation  given  by 
Eq.  (15). 


4.2.  Sensitivity  analysis 

The  following  two  cost  functions  are  considered  for  sensitivity 
analysis. 


Cost-1 :  Average  cell  voltage— Eq.  (25). 

Cost-2:  Term  responsible  for  the  concentration  polarization  at  the 
anode-electrolyte  interface— Eq.  (26). 


/i 


h 


St— Surface  area  of  the  top  surface 


In 


Ph^ 

Ph2o 


ds 


Se— Surface  area  of  the  anode-electrolyte  interface 


(25) 


(26) 


The  operating  conditions  are  the  same  as  described  in  Table  4 
and  the  applied  current  density  is  4000 Am-2.  Improving  power 
output  is  the  ultimate  goal  of  the  SOFC  design.  If  current  density 
is  fixed,  the  power  output  can  be  improved  by  increasing  the  cell 
voltage.  Sensitivity  derivatives  of  the  cost  function  representing 
the  cell  voltage  with  respect  to  various  design  parameters  can  be 
extremely  useful  in  the  design  cycle.  The  second  cost  function  rep¬ 
resents  the  term  responsible  for  the  concentration  polarization  at 
the  anode-electrolyte  interface.  As  seen,  Eq.  (26)  exhibits  an  inter¬ 
esting  nature  due  to  its  explicit  dependence  on  the  species  partial 
pressures  and  the  temperature.  Also,  improvement  in  concentra¬ 
tion  polarization  can  increase  the  SOFC  performance  and  thus,  it  is 
chosen  as  the  second  cost  function  in  the  sensitivity  studies. 

Six  design  parameters  are  included  to  compute  sensitivity 
derivatives  of  the  aforementioned  objective  functions.  The  design 
parameters  are  comprised  of  the  material  properties  of  both  the 
anode  and  the  cathode  including  porosity,  tortuosity  and  mean  pore 
radius. 

To  validate  the  implementation  of  the  discrete  adjoint  method, 
sensitivity  derivatives  computed  using  the  discrete  adjoint  method 
are  compared  with  the  same  derivatives  computed  using  the  direct- 
differentiation  method  and  the  finite-difference  method.  Note 
that  the  sensitivity  derivatives  computed  using  the  central  finite- 
difference  method  require  two  flowfield  solutions  for  each  design 
variable.  To  reduce  computational  efforts,  the  comparison  study  is 
performed  using  a  single  channel  geometry  (full  geometry  contains 
six  channels)  and  a  coarse  grid.  Relevant  physics  included  for  the 
full  geometry  is  included  in  the  comparison  study. 

Tables  5  and  6  show  the  comparisons  amongst  sensitivity  deriva¬ 
tives  computed  using  the  different  methods  for  Cost-1  and  Cost-2, 


Fig.  10.  Current  density  vectors  in  the  flowfield. 
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Table  7 

Sensitivity  derivatives  computed  using  discrete  adjoint  method. 


D.V. 

h 

h 

Sa 

-1.0110037400e-02 

-1.0899218104e-02 

Ka 

-5.0238174308e-03 

-5.1098651595e-03 

(r)a 

1.1323025342e+03 

1.1549498370e+03 

Sc 

3.3425032057e-03 

2.3809517707e-04 

Kc 

-1.6269401390e-03 

-1.2959807450e-05 

(r)c 

2.0610982993e+02 

-7.5667688361  e-01 

respectively.  As  seen,  sensitivity  derivatives  computed  using  the 
discrete  adjoint  method  and  the  direct-differentiation  method 
match  up  to  9-11  digits.  This  comparison  is  excellent.  Also,  match¬ 
ing  significant  digits  between  finite-difference  results  and  the 
discrete  adjoint  method  results  vary  between  2  and  4.  Due  to  sub¬ 
tractive  cancellation  errors,  it  is  hard  to  find  an  optimum  step  size 
when  the  finite-difference  method  is  used  to  compute  sensitivity 
derivatives.  Thus,  comparison  between  finite-difference  derivatives 
and  the  discrete  adjoint  derivatives  are  considered  satisfactory. 

After  verifying  the  consistency  of  the  derivatives  using  the 
smaller  geometry  described,  sensitivity  derivatives  have  also  been 
obtained  for  the  full  geometry,  which  includes  all  six  fuel  and  air 
channels.  Table  7  shows  the  sensitivity  derivatives  of  both  cost  func¬ 
tions  obtained  using  the  discrete  adjoint  method  for  the  original 
geometry. 

5.  Conclusions  and  future  work 

A  three-dimensional,  parallel,  unstructured  solver  has  been 
developed  to  model  complicated  transport  phenomena  present 
inside  all  components  (channels,  electrodes,  electrolyte  and  inter¬ 
connects)  of  a  solid  oxide  fuel  cell.  To  achieve  this  task,  a 
multi-species  Navier-Stokes  solver  has  been  fully  coupled  with  the 
electrochemical  solver.  The  polarization  curve  obtained  using  the 
numerical  simulation  has  been  shown  to  compare  favorably  with 
the  experimental  results  of  Wang  et  al.  [28].  Sensitivity  deriva¬ 
tives  of  two  cost  functions— cell  voltage  and  the  concentration 
polarization  have  been  computed  using  the  discrete  adjoint,  the 
direct-differentiation  and  finite-difference  methods.  Comparison 
between  sensitivity  derivatives  obtained  using  the  discrete  adjoint 
method  has  been  performed  with  that  obtained  using  other  meth¬ 
ods  to  validate  the  implementation. 

Future  work  is  targeted  at  further  developing  the  adjoint  method 
for  industrial  applications  to  fuel  cell  designs.  An  automated  design 
tool  will  be  developed  for  a  three-dimensional  SOFC  model.  A  time- 
dependent  sensitivity  analysis  will  be  implemented  to  study  the 
transient  behavior  of  the  SOFC. 
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